#陈伟毅本科毕业论文python程序1
#一维泊松方程求解

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(['science'])

# 参数设置
N = 100  # 网格点数
h = 1.0 / (N - 1)
x = np.linspace(0, 1, N)

# 定义源项和精确解
f = lambda x: np.sin(np.pi * x)
u_exact = lambda x: (1 / np.pi**2) * np.sin(np.pi * x)

# 组装刚度矩阵 A（三对角矩阵）
diag = 2 * np.ones(N) / h
off_diag = -1 * np.ones(N-1) / h
A = sp.diags([off_diag, diag, off_diag], [-1, 0, 1], format='csr')

# 处理边界条件（Dirichlet）
A[0, 0], A[0, 1] = 1, 0
A[-1, -1], A[-1, -2] = 1, 0

# 组装载荷向量 b
b = np.zeros(N)
for i in range(1, N-1):
    b[i] = f(x[i]) * h  # 简化为梯形积分近似

# 求解线性方程组
u = spla.spsolve(A, b)

# 计算误差
error_L2 = np.sqrt(h * np.sum((u - u_exact(x))**2))
error_H1 = np.sqrt(h * np.sum((np.gradient(u, h) - np.gradient(u_exact(x), h))**2))

# 可视化
fig1 = plt.figure(figsize=(4, 3))
plt.plot(x, u, 'r-', label='Numerical Solution')
plt.plot(x, u_exact(x), 'b--', label='Exact Solution')
plt.xlabel('$x$')
plt.ylabel('$u(x)$')
plt.legend()
plt.title(f'Comparison (N={N}), L2 Error={error_L2:.2e}')
plt.show()
fig1.savefig('./yiwei.png', dpi=300,  format='png')

print(f"L2误差: {error_L2:.4e}")
print(f"H1误差: {error_H1:.4e}")